%% process_showbasis
function process_showbasis(basis_type, derivative_type, varargin)
[x, y] = meshgrid(0:0.0125:1, 0:0.0125:1);
z = basis_function_ref(x(:), y(:), basis_type, derivative_type);
if isempty(varargin)
    f_begin = 1;
    f_end = size(z,2);
else
    f_begin = varargin{1};
    f_end = varargin{1};
    if f_end > size(z,2)
        f_end = size(z,2);
    end
end
triangle_region = @(x,y) (x >= 0) & (x <= 1) & (y >= 0) & (y <= 1) & (y <= -x+1);
for j = f_begin:f_end
    Z = reshape(z(:,j), size(x,1), size(x,2));
    Z(~triangle_region(x, y)) = NaN;
    surf(x, y, Z);
    hold on
end
end

%% basis_function_ref
function N = basis_function_ref(x, y, basis_type, derivative_type)
if strcmp(basis_type, "P1")
    N = zeros(size(x,1), 3);
elseif strcmp(basis_type, "P1b")
    N = zeros(size(x,1), 4);
elseif strcmp(basis_type, "P2")
    N = zeros(size(x,1), 6);
elseif strcmp(basis_type, "P3")
    N = zeros(size(x,1), 10);
elseif strcmp(basis_type, "Non-com")
    N = zeros(size(x,1), 3);
elseif strcmp(basis_type, "Hermite")
    N = zeros(size(x,1), 10);
else
    error("Invalid basis type.");
end
if strcmp(basis_type, "P1")
    if strcmp(derivative_type, "x")
        N(:,1) = 1-x-y;
        N(:,2) = x;
        N(:,3) = y;
    elseif strcmp(derivative_type, "dx")
        N(:,1) = -1;
        N(:,2) = 1;
        N(:,3) = 0;
    elseif strcmp(derivative_type, "dy")
        N(:,1) = -1;
        N(:,2) = 0;
        N(:,3) = 1;
    end
elseif strcmp(basis_type, "P1b")
    if strcmp(derivative_type, "x")
        N(:,1) = 1-x-y;
        N(:,2) = x;
        N(:,3) = y;
        N(:,4) = 27*x.*y.*N(:,1);
    elseif strcmp(derivative_type, "dx")
        N(:,1) = -1;
        N(:,2) = 1;
        N(:,3) = 0;
        N(:,4) = -27*y.*(2*x+y-1);
    elseif strcmp(derivative_type, "dy")
        N(:,1) = -1;
        N(:,2) = 0;
        N(:,3) = 1;
        N(:,4) = -27*x.*(x+2*y-1);
    end
elseif strcmp(basis_type, "P2")
    lambda = 1-x-y;
    if strcmp(derivative_type, "x")
        N(:,1) = lambda.*(2*lambda-1);
        N(:,4) = 4*x.*lambda;
        N(:,2) = x.*(2*x-1);
        N(:,5) = 4*x.*y;
        N(:,3) = y.*(2*y-1);
        N(:,6) = 4*y.*lambda;
    elseif strcmp(derivative_type, "dx")
        N(:,1) = 1-4*lambda;
        N(:,4) = 4*(lambda-x);
        N(:,2) = 4*x-1;
        N(:,5) = 4*y;
        N(:,3) = 0;
        N(:,6) = -4*y;
    elseif strcmp(derivative_type, "dy")
        N(:,1) = 1-4*lambda;
        N(:,4) = -4*x;
        N(:,2) = 0;
        N(:,5) = 4*x;
        N(:,3) = 4*y-1;
        N(:,6) = 4*(lambda-y);
    end
elseif strcmp(basis_type, "P3")
    lambda = 1-x-y;
    if strcmp(derivative_type, "x")
        N(:,1) = 0.5*lambda.*(-1+3*lambda).*(-2+3*lambda);
        N(:,4) = 4.5*lambda.*x.*(-1+3*lambda);
        N(:,7) = 4.5*lambda.*x.*(-1+3*x);
        N(:,2) = 0.5*x.*(-1+3*x).*(-2+3*x);
        N(:,5) = 4.5*x.*y.*(-1+3*x);
        N(:,8) = 4.5*x.*y.*(-1+3*y);
        N(:,3) = 0.5*y.*(-1+3*y).*(-2+3*y);
        N(:,6) = 4.5*lambda.*y.*(-1+3*y);
        N(:,9) = 4.5*lambda.*y.*(-1+3*lambda);
        N(:,10) = 27*x.*y.*lambda;
    elseif strcmp(derivative_type, "dx")
        N(:,1) = -1+9*lambda-13.5*lambda.^2;
        N(:,4) = 4.5*lambda.*(-1+3*lambda-6*x)+4.5*x;
        N(:,7) = 4.5*x.*(1+6*lambda-3*x)-4.5*lambda;
        N(:,2) = 1-9*x+13.5*x.^2;
        N(:,5) = 4.5*y.*(-1+6*x);
        N(:,8) = 4.5*y.*(-1+3*y);
        N(:,3) = 0;
        N(:,6) = -4.5*y.*(-1+3*y);
        N(:,9) = -4.5*y.*(-1+6*lambda);
        N(:,10) = 27*y.*(lambda-x);
    elseif strcmp(derivative_type, "dy")
        N(:,1) = -1+9*lambda-13.5*lambda.^2;
        N(:,4) = -4.5*x.*(-1+6*lambda);
        N(:,7) = -4.5*x.*(-1+3*x);
        N(:,2) = 0;
        N(:,5) = 4.5*x.*(-1+3*x);
        N(:,8) = 4.5*x.*(-1+6*y);
        N(:,3) = 1-9*y+13.5*y.^2;
        N(:,6) = 4.5*y.*(1+6*lambda-3*y)-4.5*lambda;
        N(:,9) = 4.5*lambda.*(-1+3*lambda-6*y)+4.5*y;
        N(:,10) = 27*x.*(lambda-y);
    end
elseif strcmp(basis_type, "Non-com")
    if strcmp(derivative_type, "x")
        N(:,1) = 2*x+2*y-1;
        N(:,2) = 1-2*x;
        N(:,3) = 1-2*y;
    elseif strcmp(derivative_type, "dx")
        N(:,1) = 2;
        N(:,2) = -2;
        N(:,3) = 0;
    elseif strcmp(derivative_type, "dy")
        N(:,1) = 2;
        N(:,2) = 0;
        N(:,3) = -2;
    end
elseif strcmp(basis_type, "Hermite")
    lambda = 1-x-y;
    if strcmp(derivative_type, "x")
        a = x.*y.*lambda;
        N(:,1) = lambda.^2.*(3-2*lambda)-7*a;
        N(:,2) = x.*lambda.^2-a;
        N(:,3) = y.*lambda.^2-a;
        N(:,4) = x.^2.*(3-2*x)-7*a;
        N(:,5) = x.^2.*(x-1)+2*a;
        N(:,6) = x.^2.*y-a;
        N(:,7) = y.^2.*(3-2*y)-7*a;
        N(:,8) = y.^2.*x-a;
        N(:,9) = y.^2.*(y-1)+2*a;
        N(:,10) = 27*a;
    elseif strcmp(derivative_type, "dx")
        b = y.*(lambda-x);
        N(:,1) = 6*lambda.*(lambda-1)-7*b;
        N(:,2) = lambda.*(lambda-2*x)-b;
        N(:,3) = -2*lambda.*y-b;
        N(:,4) = 6*x.*(1-x)-7*b;
        N(:,5) = x.*(3*x-2)+2*b;
        N(:,6) = 2*x.*y-b;
        N(:,7) = -7*b;
        N(:,8) = y.^2-b;
        N(:,9) = 2*b;
        N(:,10) = 27*b;
    elseif strcmp(derivative_type, "dy")
        c = x.*(lambda-y);
        N(:,1) = 6*lambda.*(lambda-1)-7*c;
        N(:,2) = -2*x.*lambda-c;
        N(:,3) = lambda.*(lambda-2*y)-c;
        N(:,4) = -7*c;
        N(:,5) = 2*c;
        N(:,6) = x.^2-c;
        N(:,7) = 6*y.*(1-y)-7*c;
        N(:,8) = 2*x.*y-c;
        N(:,9) = y.*(3*y-2)+2*c;
        N(:,10) = 27*c;
    end
end
end